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ABSTRACT 
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made for composite coatings of hollow cylinders. Stress 
analysis was done during a thermal cycle of successive 
heating and cooling periods. Various coating configurations 
were investigated. The results showed the effect of changes 
in coating arrangements and the failure modes due to the 
excessive stress levels. 

Solutions were obtained by approximate techniques, and 
a computer program for transient three dimensional thermal 
stress analysis was developed. The finite element technique 
was employed for the determination of the temperature 
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PART I. THERMAL STRESS ANALYSIS OF 



CYLINDERS WITH COMPOSITE COATINGS 

The efficiency of a gas turbine engine can be increased 
by raising the turbine inlet temperature. However, protec- 
tion against high temperature oxidation or sulfidation of 
the turbine blades must be provided. One method, which has 
been suggested, [1]*, is the application of ceramic coatings 
which resist high temperature corrosion. Since the material 
properties vary through such a composite blade, thermal 
stresses become an important consideration particularly 
during transients with non-uniform temperature distributions. 

The application of these ceramic coatings has been 
accomplished via a process of high-rate sputtering [1] . 

This process requires high temperatures to insure proper 
adhesion of the coatings. Since the thermal life cycle of 
the coatings differs from the thermal conditions when the 
coatings are applied, the deposition temperature must be 
considered in the design criteria. 

In order to determine the thermal stress levels within 
possible coating arrangements, hollow cylinders were analyzed 
with various material compositions. A thermal cycle con- 
sisting of a rapid heating phase followed by a cooling phase 
was used to generate temperature fluctuations and the thermal 

* Numbers in brackets refer to similarly numbered 
references in the bibliography. 
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stress distributions. The finite element technique was 
used to perform the elastic stress analysis and to determine 
the temperature distributions. 
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I. 



SIMULATION OF THE EXPERIMENTAL MODEL 



An experimental program [1] to develop the coating 
technique and to observe the coating during thermal transients 
was begun using thin walled cylinders, a simple geometry 
instead of the more complex turbine blade. Such cylinders 
were prepared by the high-rate sputter deposition of nickel 
and stabilized zirconia on a nickel superalloy substratum. 

In addition to the pure coatings, a graded composition layer, 
composed of nickel and zirconia, was used between the sub- 
stratum and the outer layer. The cylinders were exposed to 
a thermal cycle consisting of a rapid heating phase followed 
by a cooling phase. The process was accomplished by the 
cross-flow of high temperature air past the hollow cylinder, 
followed by the use of low temperature air. The analysis 
reported herein attempts to analytically simulate this 
experiment. In doing so, several simplifying assumptions 
were made. 

The gas flow around the cylinder was considered symmetric, 
therefore, only one half of the hollow cylinder was modelled. 

A plane, along the cylinder axis, passing through the forward 
stagnation point and the after stagnation point was considered 
an adiabatic surface. No axial variation of the heat trans- 
fer conditions was allowed, however, the local heat transfer 
conditions were varied in the circumferential directions. 
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Although during an actual experiment the cylinder was 
supported by some apparatus, the modelled cylinder was 
unsupported. Therefore, all transfer of heat occurred only 
in the hollow cylinder, and it was free to expand or contract 
in all directions. 

The simulation analysis was used to determine the tem- 
perature distribution and the stress distribution present 
in the cylinder at various times throughout the thermal 
cycle. Excessive stress levels were then used to predict 
when and what type of failure would occur. 

The classical thermal stress equations are coupled to 
the heat equation, however, during this simulation analysis, 
the heat equation was solved separately. The temperature 
distribution was then used to solve the thermal stress 
problem. 
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II. CYLINDER CONFIGURATIONS 



Four configurations were analyzed to determine the 
stress levels at various times in the thermal cycle. All 
geometries consisted of a hollow cylinder three inches long, 
an inside diameter of 0.5 inches, and an outside diameter 
dependent on the particular coating used. The thickness of 
the nickel superalloy substratum was 0.05 inches, and the 
thickness of the coatings varied from 0.005 to 0.009 inches. 
Drawings of all configurations are shown in Figure 1. 

A. CONFIGURATION NUMBER ONE 

The first coating investigated was a 0.005 inch layer of 
zirconia on the nickel superalloy. Material properties 
changed from those of the nickel superalloy to those of the 
zirconia coating at 0.3 inches radius. 

The finite element mesh used in this arrangement con- 
sisted of 21 isoparametric quadratic elements totaling 200 
nodes. The element discretization was specified as one ele- 
ment in the axial direction, three elements in the circum- 
ferential direction, and seven elements in the radial direc- 
tion. The seven radial elements were proportioned into four 
elements through the substratum and three elements through 
the layer. 

B. CONFIGURATION NUMBER TWO 

The second arrangement analyzed consisted of 0.002 inches 
of pure nickel on the substratum followed by 0.005 inches of 
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stabilized zirconia. The material properties abruptly 
changed at two radial positions. 

The finite element mesh for this configuration consisted 
of 21 isoparametric quadratic elements using 200 nodes. The 
subdivision in the axial and circumferential directions 
remained as in configuration number one. The seven radial 
elements were arranged with three elements in the substratum 
and two elements in each layer. 

C. CONFIGURATION NUMBER THREE 

The third arrangement analyzed consisted of 0.004 inches 
of graded layer on the substratum. This graded layer allowed 
a gradual change of the material properties from those of 
nickel to those of zirconia. This arrangement was used to 
match the material properties of the substratum to those of 
the coatings. A 0.003 inch layer of stabilized zirconia was 
placed on top of the graded layer. Properties of the graded 
composition varied in a linear fashion as described in Section 
III. 

This configuration was discretized into 24 finite elements 
totaling 226 nodes. The axial and circumferential directions 
were divided into one element and three elements respectively. 
The eight elements in the radial direction were arranged with 
four elements in the substratum, two elements in the graded 
layer, and two elements in the outer layer of zirconia. 
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D. CONFIGURATION NUMBER FOUR 



The fourth configuration was composed of the substratum 
followed by a 0.002 inch layer of pure nickel. A 0.004 inch 
coating of graded composition was placed on the pure nickel. 

A 0.003 inch coating of pure zirconia was used as the outer 
layer. 

The finite element mesh for this configuration consisted 
of 27 isoparametric quadratic elements using 252 nodes. The 
axial and circumferential directions were subdivided as in . 
the other configurations. Nine elements were used in the 
radial direction, of which three were placed in ' the substratum 
followed by two in each of the three layers . 

E. COMMENTS ON THE FINITE ELEMENT MESH 

Emphasis has been placed on the radial and tangential 
variations, therefore, the mesh was made finer for those 
directions. Since a limitation on the computer core require- 
ment and computing time dictates the number of elements and 
nodes, the axial direction was modelled using only one element 
The problem solution is a three dimensional method, there- 
fore, ideally it would be more appropriate to have additional 
elements. Since the main interest was focused on the radial 
variation in material composition and the circumferential 
variation of heat transfer coefficients, the mesh used was 
considered appropriate. There was no axial variation in 
material properties or boundary conditions, therefore, no 
axial variation in the temperature distribution resulted. 
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A slight variation in the stress distribution resulted due 
to the three dimensional nature of the problem. 
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FIGURE 1. CYLINDER CONFIGURATIONS 
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III. MATERIAL PROPERTIES 



The materials used in the analysis were a nickel-super- 
alloy (IN-738) , pure nickel, and stabilized zirconia (Zr O 2 
plus Ca 0) . The substratum of the hollow cylinder was com- 
posed of the nickel superalloy, and the coatings were com- 
binations of nickel and zirconia. Material property values 
were assumed constant throughout the temperature range inves- 
tigated. The values used were the average of those available 
in the literature [2,3]. The numerical values for elastic 
modulus , Poisson ' s ratio, coefficient of linear expansion, 
thermal conductivity, and volumetric heat capacity (density 
times specific heat) of the various materials are given in 
Table 1. 

Matching of material properties was attempted by using 
a graded composition between layers in some of the coating 
configurations. Approximation to the properties of these 
compositions was made by using a function based on the 
radial position. If two materials A and B have properties 
P^ and Pg, the property of the composition P c was determined 
by 



r B " r r " r A 

p = p (_ ) + p ( _ 

C A v r^ - r, r B r_ - r. 



’B 



'B 



(1) 



where r A is the radius of material A and r B , that of material 
B. 



19 



TABLE 1 



MATERIAL PROPERTIES 



Nickel Superalloy - IN-738 

Modulus of Elasticity - 

Poisson's Ratio 

Coefficient of 
Linear Expansion 

Thermal Conductivity 

Volumetric Heat Capacity 

Pure Nickel 

Modulus of Elasticity 

Poisson's Ratio 

Coefficient of 
Linear Expansion 

Thermal Conductivity 

Volumetric Heat Capacity - 

Stabilized Zirconia 

Modulus of Elasticity 

Poisson's Ratio 

Coefficient of 
Linear Expansion 

Thermal Conductivity - 

Volumetric Heat Capacity - 



31.75 x 10 6 PSI 
0.3 

8.6 x 10 ^ in/in-deg.F 
1.34 x 10 ^ BTU/min-in-deg.F 
0.029 BTU/cu.in.-deg.F 

30.00 x 10 6 PSI 
0.3 

8.6 x 10 ^ in/in-deg.F 
6.89 x 10 ^ BTU/min-in-deg.F 
0.035 BTU/cu.in.-deg.F 

24.8 x 10 6 PSI 
0.3 

5.0 x 10 ® in/in-deg.F 

3.36 x 10 ^ BTU/min-in-deg.F 
0.032 BTU/cu.in.-deg.F 
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The materials used in the graded layer were pure nickel 
and stabilized zirconia. The values for the various pro- 
perties of the graded compositions are shown in the material 
property graphs, Figure 2. The property values used during 
the quadratic variations were determined by specifying 
property values for three points in the layer and then forming 
the quadratic function. Two of the points were specified 
to have properties of the inside and outside of the graded 
layer, i.e. , the properties of nickel and zirconia. The 
midpoint was specified as the arithmetic mean of the pro- 
perties at the midpoint of the linear variation and the 
properties of zirconia. The resultant quadratic variation 
is shown on the same plots as the linear variation. 
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(*) Configuration Three: a=0.3, b=0.302, c=0.304 

Configuration Four : a=0 . 302 , b=0.304, c=0.306 
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FIGURE 2. MATERIAL PROPERTY GRAPHS 
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IV. THERMAL CYCLE CONDITIONS 



The conditions which specify the thermal stress problem 
are classified as boundary conditions and initial conditions. 
With respect to the temperature distribution problem, the 
boundary conditions consist of convection heat transfer 
coefficients and the external fluid temperatures, and the 
initial condition is the initial temperature distribution. 

The stress problem has boundary conditions which specify 
zero or given displacement at various positions. In addi- 
tion, the zero-stress temperature distribution must be 
supplied since all temperature distributions are compared 
to the zero stress condition. 

A. HEAT TRANSFER BOUNDARY CONDITIONS 

The thermal cycle consisted of a rapid heating phase and 
a cooling phase. The cylinder wall was heated from room 
temperature (77 deg. F) to approximately 2000 deg. F in a 
period of 45 seconds. During the cooling phase, the heat 
transfer conditions were changed to cool the cylinder back 
to room temperature in approximately two minutes. During 
the heating phase, a convection boundary condition was main- 
tained on the external surface and an adiabatic condition on 
the internal surface. The cooling phase was characterized 
by convection boundary conditions on both internal and 
external surfaces. A circumferential variation of the 
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external heat transfer coefficients was used during both 
phases. The internal heat transfer coefficient was uniform 
during the cooling phase. No axial variation in the heat 
transfer coefficients was used. 

In order to determine an average heat transfer coeffi- 
cient for both the heating and cooling phases of the transi- 
ent cycle, the cylinder was treated thermally as a lumped 
capacitor. Using this average coefficient, an average 
Nusselt number was calculated using air as the external 
heating fluid. The temperature of the external fluid was 
assumed to be 2500 °F. A Reynolds number for flow past a 
cylinder was then computed using the Hilpert correlation. 

This Reynolds number was used to determine the local Nusselt 
number and therefore heat transfer coefficient for air flow 
past a horizontal cylinder. These values of Nusselt number 
were obtained from Kreith [5] and Krall and Eckert [6]. 

For the cooling phase, the average heat transfer coefficient 
was used for the internal surface. The values used during 
the computer analysis are shown in Table 2. These values 
were used for all configurations. 

The resultant temperatures of the cylinder wall were 
considered adequate for the analysis. The temperature 
variation at 45 seconds was between 1850 deg.F and 2050 deg F. 
Room temperature (77 deg. F) was reached after 2.25 minutes 
into the cooling phase. 
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TABLE 2 



Heating Phase 



Fluid Temperature: 2500 deg. F 

Average Heat Transfer Coefficient: 22.6 BTU/hr-f t -deg. F 

Average Nusselt Number: 22.5 

Reynolds Number: 2262 

Equivalent Fluid Velocity: 124.6 ft/sec 

Local Heat Transfer Coefficients 

Angle Position from Stagnation Coefficient 

(degrees) (BTU/hr-f t2-deg . 



0.0 


40.0 


30.0 


41. 0 


60.0 


38.0 


90.0 


26.0 


120.0 


' 11.5 


150.0 


16.0 


180.0 


12.0 



Cooling Phase 



Fluid Temperature: 60 deg. F ~ 

Average Heat Transfer Coefficient: 13.5 BTU/hr-f t -deg. F 

Average Nusselt Number: 47.1 

Reynolds Number: 7300 

Equivalent Fluid Velocity: 22.9 ft/sec 



Local Heat Transfer Coefficients 



Angle Position from Stagnation 


Coefficient 


(degrees) 


(BTU/hr-f t2-deg. 


0.0 


22.0 


30.0 


20.5 


60.0 


15.0 


90.0 


9.0 


120.0 


6.0 


150.0 


12.0 


180.0 


12.0 



F) 



F) 
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B. STRESS ANALYSIS BOUNDARY CONDITIONS 



The hollow cylinders analyzed were free to expand or 
contract in all directions. 

C. INITIAL CONDITIONS 

The cylinder was specified to be at an initial uniform 
temperature of 77 deg. F. The zero-stress temperature 
distribution was assumed to be the temperature when the 
coatings were applied. This was specified to be 932.0 deg. F. 
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V. RESULTS OF THE ANALYSIS 



The results of the analysis showed that the stresses 
within the surface coatings varied from compression at room 
temperature to tension at the maximum temperature. Due to 
the variation of the heat transfer boundary conditions, a 
uniform temperature equal to the zero-stress condition was 
not present during the thermal cycle. The complex nature 
of the circumferential variation in heat transfer coeffi- 
cients causes temperature and stress distributions vastly 
different from those common to radial heat flow. 

All configurations follow the same general stress history. 
During the period when the temperatures are below the zero- 
stress condition, the coatings remain in compression. At 
the same time, the sign of the substratum stresses varies 
according to the position. At approximately 12 seconds, the 
temperature distributions reached the zero-stress condition 
and the stresses in the coatings became tensile. These ten- 
sile stresses increased as the temperature increased 
reaching a maximum at 0.75 minutes. All configurations 
followed this general scheme. The maximum tensile stress 
is dependent on the zero-stress condition, i.e., the higher 
the zero-stress temperature, the lower the maximum tensile 
stress. Configuration Four had somewhat lower stress levels 
than the other configurations; however, the difference is 
not considered significant. 
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The use of a graded layer was investigated further by 
specifying a quadratic composition instead of a linear 
composition. An adjustment of the thermal conductivity and 
the coefficient of linear expansion to follow a quadratic 
function produced only minor changes in the stress levels. 
When all properties were adjusted to the quadratic composi- 
tion variation, no appreciable change occurred in the 
stress levels. 

Although the maximum stress levels are excessive, the 

-3 

maximum strains (9.6 x 10 ) are considered small enough to 

justify the use of a linear elastic analysis. In addition, 
the brittle nature of zirconia justifies the use of a brittle 
fracture failure criteria. Therefore, the elastic analysis 
provides a reliable method for determining whether the struc- 
ture will fail or not. Using the maximum principal stress 
theory of failure, some prediction of the time and tempera- 
ture when failure occurs can be made. If the ultimate 
strength of the ceramic material is set at 18.0 KSI and 
that of the nickel superalloy at 90.0 KSI, observation of 
the average tangential and axial stress levels for all con- 
figurations indicates that failure would occur at about 0.2 
minutes and an average temperature of 1000.0 °F. At most, 
this is only a rough calculation. Observations of an actual 
specimen could confirm this estimate by raising and lowering 
the specimen to successively higher temperatures, followed 
by appropriate microscopic examination. 
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The graphical representations of the thermal stress 
and temperature distributions were made using the numerical 
values obtained for the nodes at the cylinder ends. No 
attempt was made to smooth out the curves and no interpola- 
tions were performed. Each plot is identified by a time and 
a configuration designation. For example, the first plot is 
identified by TIME = 0.05 CONF.l-O. This plot represents 
the conditions at 0.05 minutes after start of the thermal 
cycle for configuration number one at position zero degrees 
relative to the stagnation point. Each axis has been labelled 
and the units identified. Since the range of values for 
stress and temperature varies from one graph to another, the 
scale factors used in plotting are provided on each plot. 

In order to provide an overall view of the stress levels 
throughout the thermal cycle, a plot of extreme principal 
stresses versus time has been prepared for each configuration, 
Figures 7, 14, 21, 26. 

A. CONFIGURATION NUMBER ONE 

This configuration contained a single layer of stabilized 
zirconia on the nickel superalloy substratum. Comparison 
of the maximum principal stresses in all four configurations 
shows that configuration number one has the highest stress 
level. This maximum stress level (140.0 KSI) occurred at 
the maximum temperature, i.e., at the end of the heating 
phase, and was located at a radius of 0.3017 inches, the 
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zero angle position, and at the midpoint of the cylinder 
length. The largest tensile stresses present in configura- 
tion number one were all tangential and axial components 
of the stress state. The following figures show the radial 
variation of temperature and stresses at a given angular 
position for selected times during the heat transfer cycle: 
Figure 3 (a through j) at zero degrees, Figure 4 (a through j) 
at 60 degrees, Figure 5 (a through j) at 120 degrees. Figure 6 
(a through j) at 180 degrees. (Similar sets of figures are 
provided for the remaining three configurations.) 
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B. CONFIGURATION NUMBER TWO 



This configuration contained a layer of pure nickel 
between the substratum and the layer of zirconia. The 
maximum principal stress (134.9 KSI) occurred at the maximum 
temperature condition and was located at a radius of 0.3045 
inches, the zero angle position, and at the midpoint of the 
cylinder length. As with the other configurations, the 
larger tensile stresses were in the tangential and axial 
directions . 

The comparison of stress and temperature distribution 
as a function of angle position is shown in Figure 12. The 
axial stress level is nearly constant as the circumferential 
position changes, however, the radial and circumferential 
stresses change at about 90 - 120 degrees. This change in 
the stress level coincides with a change in the temperature 
distribution. The change in the temperature distribution 
results from the changes in the temperature boundary conditions 
in the down stream portion of the air flow past the cylinder. 
Figure 12 was prepared using data for the cylinder ends at a 
radius of 0.3045 inches and a time of 0.75 minutes. 

The time dependence of the stress levels is shown in 
Figure 13. The tangential and axial stresses for a position 
on the cylinder end at a radius of 0.3045 inches and an angle 
of 60° vary in the same manner as the temperature at that 
point. The radial stress is minimal throughout the time 
period. 
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C. CONFIGURATION NUMBER THREE 



This arrangement of coatings included a graded layer of 
nickel and zirconia. The effect of this layer was to 
eliminate the abrupt change in the stress levels seen in the 
first two coating arrangements. The maximum principal stress 
(134.9 KSI) occurred at the maximum temperature. The minimum 
principal stress (-97.2 KSI) occurred at the beginning and 
end of the thermal cycle, i.e., the room temperature condition. 
Since the tangential and axial stress levels are excessive 
for a major part of the cycle, failure would occur due to 
these stresses. 

As with configuration two, a comparison of stress and 
temperature distributions as a function of angular position 
was prepared using data for the cylinder ends at a radius 
of 0.3055 inches and at a time of 0.75 minutes. The graphical 
representation of this comparison is shown in Figure 19. The 
function of stress and temperature versus time for the 60° 
angular position, at a radius of 0.3005 inches at the cylinder 
ends is shown in Figure 20. 
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FIGURE 15(a) 
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FIGURE 15(b) 
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D. CONFIGURATION NUMBER FOUR 



This coating arrangement included a layer of pure nickel 
between the substratum and the graded composition of 
configuration number three. The maximum principal stress 
(130.3 KSI) which occurred at the maximum temperature 
condition, was located at a radius of 0.3075 inches, the 
zero angle position, and at the midpoint of the cylinder 
length. As with the other configurations, the tangential 
and axial stresses were the largest for this configuration. 
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FIGURE 25 (i) 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



The thermal stress analysis of the various cylinder 
configurations showed that the maximum stress levels of the 
four coating arrangements were approximately equal. The 
maximum stress level was excessive enough to cause failure 
of the coatings and was dependent on the deposition temp- 
erature and the maximum temperature distribution. The abrupt 
change in the stress levels between the substratum and the 
outer coatings of configurations one and two was reduced by 
the use of the graded layers in configurations three and 
four. 

The following recommendations are suggested for future 
development and analysis. 

A. MATERIAL PROPERTIES 

An extensive investigation of the temperature dependent 
properties should be completed. In addition to the various 
individual materials, the graded layer properties must be 
thoroughly understood before a complete design analysis 
can be completed. All creep data for these materials must 
be obtained. If possible, material properties for the 
compressive stress state should be obtained. 
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B. VARIATIONS TO THE PROBLEM 



The following situations should be investigated. 

1. Other geometries, such as wedges and turbine blades, 
would exhibit different temperature and stress distributions. 
Investigation of these geometries may show lower stress 
levels which would be less susceptible to failure. These 
investigations should be made. 

2. The temperature distribution and maximum stress levels 
are directly related to the thickness of the zirconia coating. 
This coating depth should be varied to determine an optimal 
thickness for the minimization of the stress levels. 

3. Since the deposition temperature affects the maximum 
tensile stresses, more investigation to determine the most 
appropriate temperature for coating application should be 
made. More information about the compressive state will 

also allow a better evaluation of the effect of the deposition 
temperature selection. 

4. Although changes in the material property variation 
of the graded layer produced minor changes in the stress 
levels, the effect of the thickness of the graded layer should 
be investigated. 

5. Since thermal load conditions are only part of the 
total load conditions on a turbine blade, an investigation 

of the effect of body forces on the stress levels in the blade 
should be completed. The program used in this analysis allows 
the specification of these transient body forces. 
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PART II. TRANSIENT THERMAL STRESS ANALYSIS 



The classical equations governing thermal stress analy- 
sis are coupled to the heat conduction equations, however, 
in most studies of this topic including this one, the stress 
analysis and the heat conduction are separated. The analysis 
of the transient thermal stress condition is performed in 
two steps. In the first one, the unsteady heat equation is 
solved for the temperature distribution. Then, using this 
temperature distribution, the equations of elasticity can be 
solved to determine the stress distribution. Analytically 
this solution process is possible only in the simplest of 
problems. However, through the use of the finite element 
techniques, an approximate solution can be obtained for more 
realistic problems. 

The computer program which was used as a basis to solve 
the heat conduction equation was written by Lew [7] . Several 
modifications were made to this program as described in 
Section I. The computer program used as a basis to solve 
the stress problem was written by Leonidas [8] . Several 
changes, as described in Section II, were also made to this 
program. 

Since three dimensional finite element analysis generates 
such a voluminous amount of information, an additional post- 
processing program was written to condense the relevant 
information . 
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All programs were written using Fortran IV 
language and solutions were obtained using the 
computer system. 



x 



G Level 
IBM 360/67 
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I. 



TEMPERATURE DISTRIBUTION 



The unsteady heat conduction equation, expressed in 
rectangular coordinates, is 



£» + w (k y ff> + ^ <k z fi> + ° - » c fi 



3T. 



3T 



(2) 



This equation is subject to the following boundary conditions 



T = T_. on the boundary or 

r> 



k |2 i + k |I 

x 3x x y 3y 



1 + 

y 



3T 
3 z 



1 + 
z 



+ h (T-T f ) = 0 



on the boundary 



(3) 



and a specified initial temperature distribution. The 
symbols used in these equations have the following definitions 



Q 






- rate of internal energy generation per 
unit volume 

- density 

- specific heat 

- thermal conductivity in the x,y, and z 
directions 

- the direction cosines at a point on the 
boundary 

- external heat flux 

- surface heat transfer coefficient 

- fluid temperature 
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As described in Ref. 9, the unsteady heat equation can 
be transformed into a discretized form by the finite element 
technique. Using the symbol [ ] for an n x n matrix, 

< > for a row vector, and { } for a column vector, 

[H] {T} + [C] {T} + {F} = 0 (4) 

where the temperature T of equation (2) has been expressed 
as the product of the nodal shape functions and the nodal 
temperature values . 



T = <N> {T} 



(5) 



3T 

The time derivative of the temperature — has been expressed 

dt 

as 



8T 

3t 



T 



<N> 



{T} 



(6) 



The matrices [H] , [C] , {F} are defined at the element 

level by 



!H) e 

tC1 e 



V 



,<k {M)<M> 

J x L 3x 3x 



y 3y < 3y > 



+ k 



z { H }< H >)dx *y 



/{N} pC <N> dx dy dz 
V 

e 



dz 

(7) 

(8) 
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{F} = - / Q {N }dx dy dz + /q{N}dS + ( / (N}h<N>dS ) {T. } 

6 V S e 1 

( 9 ) 



where 



q = q - h T, 



The program developed by Lew [7] solves the discretized 
equation subject to the listed boundary conditions. In order 
to adapt the program to the thermal stress problem, several 
alterations were incorporated. 

A. TIME INTEGRATION SCHEME 

Reference 7 provides numerous time integration schemes 
to deal with the first time derivative of the temperature. 

In this study only the best integration schemes were retained. 
A five point backward difference formula was used to approxi- 
mate the time derivative. For the ith time step, this 
approximation can be expressed as 



(|?) • = T. = r4-(3T. . - 16T . , + 36T. _ - 48T. + 25T. ) 

3t i i 12At i-4 1-3 i-2 l-l i 



( 10 ) 



where At is the specified time step. The truncation error 

1 4 3^T 

for this expression is +r(At) — 

at 

If equation (10) is written in vector form and substituted 
into equation (4) , the final form of the discretized heat 
equation is obtained. 
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(12At [H] + 25 [C] ) • {TV} = 

48{T i-l } “ 36{T i-2 } + 16{T i-3 } " 3{T i-4 } “ 12 A t{F } 

(ID 

If [H] , [C] , {F}, At and the four preceding values of {T\ } 

are known, the system of simultaneous equations is solved 
for the next set of {TV }. 

The initial four values for the temperature distribution 
is obtained by using the trapezoid method for a specified 
number of steps at a small At. After each n steps at At, 
the values of {T^ } are saved. When three sets of {T^} plus 
the initial known distribution are available, the transient 
solution using a time step of (nAt) can proceed. For exam- 
ple, an initial step size of 0.001 taken for 150 steps 
results in a step size of 0.05 for the remainder of the 
transient problem. 

Using the solution technique described above, the trans- 
ient temperature distribution for an infinite hollow cylinder 
were determined. The problem characteristics are: 



Inside diameter 

Outside diameter 

Initial temperature 

Internal and external 
heat transfer coefficients 

Internal Fluid temperature 

External Fluid temperature 

Initial time step size 

Number of initial steps 



= 0.5 inches 

= 0.6 inches 

= 932 deg. F 

= 25.0 BTU/hr-sq. ft .-deg .F 
= 896 deg. F 

= 1994 deg. F 
= 0.001 min. 

= 150 
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Results of the approximate solution and the exact solution 
are given in Table 3. 

B. NODAL MATERIAL PROPERTIES 

The usual method of specifying properties to be constant 
throughout a finite element becomes troublesome when com- 
posite materials have continuously varying properties. Many 
separate elements would be required in order to adequately 
represent a structure which has material properties that are 
continuous functions of spatial coordinates. 

In order to alleviate the problem, several routines in 
the basic program of Ref. 7 were altered to accommodate 
nodal properties in addition to element properties. A 
property P is expressed as P = N. P., where P^ is a vector 
of nodal values and N^ are the same shape functions used in 
equation (5) for the temperature variable. The matrices 
[H] , [C] , and {F} are formed by the same integration tech- 

niques as used in the basic program. However, the proper- 
ties are no longer constant during the integration, and 
must be re-evaluated at each integration point. 

As an example, consider steady state heat flow through 
a cylinder in which the thermal conductivity varies contin- 
uously along a radius. As proposed above, the development 
of a twenty nodal point element could represent exactly any 
linear or quadratic variation of material properties. Table 
4 displays the results of such a problem as compared to 
the exact solution. The column heading "Element Property 
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TABLE 3 



TRANSIENT TEMPERATURE OF THE INFINITE HOLLOW CYLINDER 

Exact Solution (Deg. F) 



Radius 

(in) 


t = 0.25 
(min) 


0.5 


1.0 


3.0 


0.25 


1277.165 


1409.796 


1479.776 


1491.715 


0.26 


1278.061 


1410.916 


1481.013 


1492.977 


0.27 


1279.104 


1412.062 


1482.214 


1494.182 


0.28 


1280.291 


1413.236 


1483.381 


1495.348 


0.29 


1281.618 


1414.438 


1484.518 


1496.473 


0.30 


1283.081 


1415.669 


1485.626 


1497.561 



Finite Element Solution 





Using 5-Point 


Finite Difference Method (Deg. 


F) 


Radius 


t = 0.25 


0.5 


1.0 


3.0 


(in) 


(min) 








0.25 


1277.220 


1409.897 


1479.879 


1491.810 


0.26 


1278.116 


1411.016 


1481.116 


1493.068 


0.27 


1279.160 


1412.163 


1482.316 


1494.278 


0.28 


1280.347 


1413.336 


1483.483 


1495.443 


0.29 


1281.674 


1414.538 


1484.612 


1496.568 


0.30 


1283.138 


1415.769 


1485.726 


1497.654 
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TABLE 4 

STEADY HEAT FLOW THROUGH A HOLLOW CYLINDER 



Inside 


Diameter 


0.5 in. 




Outside 


Diameter 


0.6 in . 




Inside 


Surface Temperature 


500.0 °F 




Outside 


Surface Temperature 1000.0 °F 




Thermal 


Conductivity Variation in Outward Radial 


Direction : 


(a) 


For linear function 


, 0.05 to 0.25 BTU/in- 


min-°F 


(b) 


For quadratic function, 0.2025 to 0.225625 




(at mid-radius) to 


0.25 BTU/in-min-°F 




Linear Variation of Thermal Conductivity 






(Two Radial 


Elements) 






Exact 


Element Property 


Nodal Property 


Radius 


Solution 


Solution 


Solution 


(in.) 


(deg. F; 


(deg. F) 


(deg. F) 


0.25 


500.0 


500.0 


500.0 


0.2625 


575.884 


590.871 


576.502 


0.275 


669.145 


677.510 


670.301 


0.2875 


794.697 


842.256 


796.982 


0.3 


1000.0 


1000.0 


1000.0 


Quadratic Variation of 


Thermal Conductivity 






(Two Radial 


Elements) 






Exact 


Element Property 


Nodal Property 


Radius 


Solution 


Solution 


Solution 


(in.) 


(deg. F) 


(deg. F) 


(deg. F) 


0.25 


500.0 


500.0 


500.0 


0.2625 


623.879 


627.602 


624.213 


0.275 


748.213 


749.253 


748.803 


0.2875 


873.441 


877.342 


873.747 


0.3 


1000.0 


1000.0 


1000.0 
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Is? 



Solution" identifies the results obtained by specifying 
properties constant over an element. The column heading 
"Nodal Property Solution" identifies the results obtained 
by specifying the properties variable over an element, 
in other words nodal properties . 

C. BOUNDARY CONDITIONS 

The load vector represented by equation (9) contains 
two surface integrals which pertain to boundary conditions 
described by equation (3) . The development of this load 
vector is shown in Ref. 9. 

The usual technique for specification of the boundary 
conditions is to consider q and h as constants over the 
applicable element surface. In order to specify convection 
boundary conditions which vary over an element surface, the 
flux and heat transfer coefficient have been treated as 
nodal values in the same manner as nodal material properties . 
The heat generation term, Q, is also treated in this manner. 
Therefore, instead of the constant terms in equation (9) , the 
usual shape functions N^ are used with the nodal values to 
form the variable terms Q = N^ Q^, q = N^ q^, and h = N^ h^ . 
The integration processes proceed as described in Ref. 7, 
with the addition of these terms. 

The surface integrals are computed using a coordinate 
transformation to allow the use of Gauss Quadrature. This 
transformation maps the Cartesian coordinates (x,y,z) into 
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a local coordinate system (£,ri/C). As described by Hanson 
[ 10 ], the unit area vector in the (x,y,z) system is evaluated 
in the (£,n,C) system by 



" ds e = * |f> d “ « , ( 12 ) 



where 



r = x(£,rwc)i + y(5/n»c)j + z(£,ri/C) k (13) 



specifies the surface and a, 3 represent £,ri/ or 

3 r 3 r 

The vector V = -r— x -r-r- is normal to the surface area 

da 83 

da dg and must be evaluated during the integration process. 
It is convenient to express equation (12) as 



n dS = | V | (- 3 -) da dB = |V|n da dB , (14) 

l v l 

A 

where n is a unit normal vector to surface area da d$ . 
Therefore, during the integration process, the factor 




t<|* 

8a 



82 

86 



ll IX) 

8 a 86 



, 8 x 8 z 
8 a 86 



8 z 8 x , 
8 a 86 ; 



+ 



. * * „ 2 1/2 
/ 8x 8y _ d£ 3x> 1 

' 8a 86 8a 88 



(15) 



must be evaluated. If the surface integral is performed 
over an area on which one of the variables x, y, or z is 
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constant, the quantity |v| becomes identical to the deter- 
minant of the Jacobian matrix for the surface transformation. 
For a general surface r(x,y,z), the complete expressions for 
|V| must be used. 



II. STRESS DISTRIBUTION 



The finite element technique provides the required 
methods for the solution of the three dimensional elastic 
stress analysis problem. Reference 9 describes the deri- 
vations for the following equations. A structure can be 
divided into several elements, each formed by a specified 
number of nodes. If this structure is acted upon by various 
forces, the nodal displacements {r} are determined by 
solving the following equation: 

[K] {r } = {R} - {F } - { F}_ - {F} - {F} , (16) 

B P £ o 

where the matrices and vectors are described by 



[K] 


= / 

V 


[B] T 


[DJ [B] dV 


(17) 


[f) b 


= -/ 

V 


[n] t 


{ V dv 


(18) 


[F] p 


= -/ 

s 


[N) T 


{f } dV 
P 


(19) 


IFlc 

o 


= -/ 

V 


[B) T 


[D] {e o > dV 


(20) 


*■3 

Q 

0 


= +/ 

V 


[B) T 


{a } dV 
o 


(21) 
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and { R} is the vector of concentrated nodal loads. The 
other symbols used in these equations are defined as: 



[N] 


- vector 


of shape functions 


[B] 


matrix 


of shape function derivatives 


[D] 


elasticity matrix 


< f B> 


vector 


of body forces 


<V 


vector 


of surface forces 


< e o» 


- vector 


of initial strains 


<«o> 


vector 


of initial stresses 



Following the solution of equation (16) , the nodal displace- 
ments {r} are used to evaluate the nodal strains. 

{ e } = [B] {r} (22) 

The final step in the analysis is the calculation of the 
nodal stresses. 



{ 0 } = [D] ({ e } - { £ Q } ) + (a 0 ) (23) 

A computer program to make these calculations was 
written at the Naval Postgraduate School by Leonidas [8]. 

In that version of the program, the right hand side load 
vector was allowed to contain only concentrated loads (R). 

In order to permit the presence of thermal stresses and body 
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forces, several changes were made to allow {e } and {f D } 

O ID 

to have non-zero values. 

A. GENERATION OF LOAD VECTORS 

Temperature differences always generate strains but 
sometimes stresses may still be zero, for example, when 
the temperature differences are constant throughout, e in 
equation (23) would be equal to giving zero stress for 
one material. The presence of a non-uniform temperature 
distribution in a solid body causes stresses and strains due 
to non-uniform expansion. In addition, a solid body composed 
of materials with differing coefficients of linear expansion, 

A, will at times experience stresses and strains. This happens 
when the temperature distributions differ from that of the 
zero-stress condition. 

The forces present in a body during a thermal stress 
state can be calculated by setting = AAT, where AT is 
the difference between the actual temperature distribution 
and the zero-stress temperature distribution. In order to 
form a consistent load vector, these forces, calculated by 
equation (20) , are added to any other nodal forces present 
on the right hand side of equation (16) . 

The body forces, due to some force field, are converted 
into nodal forces by using equation (17) . These are added 
to the right hand vector of equation (16) to form the total 
consistent load vector. 
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B. NODAL MATERIAL PROPERTIES 



As discussed in Section IB, the use of nodal material 
properties is a more convenient method for composite body 
analysis. In the stress analysis operations, the elastic 
matrix and the coefficient of linear expansion enter into 
the calculations. The elastic matrix is composed of the 
elastic modulus and Poisson's ratio. For the different 
operations, these properties are evaluated by using 
E = N^ E^, y = N^ vu , and A = N^ A^ . During the various 
integration processes these property values are evaluated 
at each integration step. 

The use of nodal properties requires careful considera- 
tion as to the effect the shape functions have on the pro- 
perties within the element. If an element of one material 
is adjacent to an element of another material, the use of 
nodal properties can lead to erroneous results. Since the 
nodes common to the elements can have only one set of pro- 
perty values, one of the elements will not have the expected 
uniform properties. It would be better to use constant ele- 
ment properties or combine the two elements into one and 
use variable properties over that element. The amount of 
property variation between the two elements will also assist 
the analyst in determining the proper technique. 
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III. DESCRIPTION OF THE PROGRAMS 



Three programs were used to complete the thermal stress 
analysis. As noted in Sections I and II, the programs of 
Lew and Leonidas [7,8] were used as a basis to develop the 
programs for the temperature distribution and elastic stress 
analysis. The third program was written to process the 
information generated in the stress analysis. The computer 
listings for these three programs are provided in Appendices 
F, G, and H. 

A. TEMPERATURE PROGRAM 

The main alterations to the temperature program [7] were: 
(1) the elimination of the second time derivative of the 
quasi-harmonic equation, (2) the addition of the time inte- 
gration scheme described in Section IA, and (3) changes to 
subroutines FORMC , CUBE, and BCOND to allow the specification 
of nodal material properties and nodal boundary conditions. 
The version of the program listed in Appendix F allows 278 
nodes and a bandwidth of 40. These specifications were 
varied to adjust the computed core requirements to the 
particular problem being solved. A problem with 278 nodes 
and a bandwidth of 40 required approximately 350,000 bytes 
of storage. 

The assumptions made in the solution of a problem using 
this program are: (1) material properties are constant with 

time and independent of temperature, (2) directional material 



* 
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properties must align with the global coordinate axes, and 
(3) the boundary conditions are constant within a heat 
transfer phase, but may be altered at the end. 

1 . Subroutine TIMEHB 

As discussed in Section IA, this subroutine performs 
a time integration using a five point backward difference 
method. During each heat transfer phase, boundary conditions 
are constant. Following applications of these boundary con- 
ditions to the appropriate matrices and vectors, the trape- 
zoid integration method is used to generate three starting 
values. Along with the initial temperature distribution, 
these three starting values provide the necessary information 
to use the five point difference scheme. The output of tem- 
perature distributions is controlled in a manner similar to 
that used in the original program [7] . These output values 
are adjusted to the reference temperature supplied in the 
data deck. During a stress analysis problem, this reference 
temperature is the zero-stress temperature. Since these 
temperature distributions are used in a subsequent stress 
analysis, proper control is provided to punch cards or 
transfer the information to temporary storage devices. The 
program in its present form allows up to twenty separate 
temperature vectors to be stored for future use. 

2. Subroutine FORMC, CUBE 

Subroutine FORMC was modified to allow nodal proper- 
ties to be used in the integration process for the matrices 
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[H] and [ C ] . The subroutine CUBE performs the cubic inte- 
gration used in the formation of these matrices and there- 
fore required slight modification. 

3 . Subroutine BCOND 

Subroutine BCOND and its subsidiary routines were 
modified to accommodate nodal boundary conditions. Subrou- 
tine JACOB was altered to evaluate the factor given in 
equation (15). This factor was used in the surface integra- 
tion as discussed in Section IC. 

4 . Data Preparation for the Temperature Program 

The data required for a transient solution is com- 
prised of eight groupings of cards. Careful detail must be 
followed to insure proper formats are used and the correct 
order of cards maintained. A data check can be performed 
without any major calculations. The data check normally 
requires 20-30 seconds, however, the full computer core 
requirement is needed. 

Preparation of the data cards is explained in the 
initial comment cards of the program listing in Appendix 
F. A sample data deck for a 32 node, 2 element program 
is provided in Appendix A. Several aspects of the data 
preparation must be thoroughly understood. 

a. Nodal Properties 

The method of input for values which are common 
to several nodes was provided to reduce the number of input 
cards. Care must be taken to insure that every node has 
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properties assigned to it, and that the order of material 
property cards follow the order of ascending nodal numbers. 

b. Output Storage 

The number of times stored output is requested 
must match the number of loads used in the stress analysis 
program. 

c. Steady State Solution 

If a steady state problem is being analyzed and 
no convection or internal generation conditions exist, a 
blank card must be inserted in place of the heat transfer 
phase data. 

d. Stacking of Problems 

If this program is used in conjunction with the 
elastic stress analysis program, only one problem can be 
loaded at a time. Therefore, the final comments concerning 
stacking problems do not apply, and a "STOP" card must be 
used. 

e. Boundary Condition Data 

When preparing the boundary condition data, the 
number of cards with boundary property values (NODECD) is 
given on the first card of each boundary condition data 
set. Ensure that NODECD cards are used in that particular 
set. 

5 . Job Control Language 

The necessary job control language for using the 
temperature program is listed in Appendix B. The two 
cards beginning with //GO.FTlOFOOl and //GO . FT11F00 1 are 
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necessary for the storage of the matrices [H] and [C] during 
the time integration process. They are not required for 
the steady state solution. The card beginning with 
//GO . FTO 7F0 01 is used to store temperature distributions 
on temporary files for use in the stress analysis program. 

If punched output is desired, this card should be replaced 
by //GO. FTO 7 F0 01 DD SYSOUT=B 

The card beginning with //GO.FT06F001 is necessary to allow 
large amounts of output from the computer printing devices . 

If less than 3000 lines of output is anticipated, this card 
can be removed. These job control cards are designed for 
use at the W.R. CHURCH COMPUTER CENTER, NAVAL POSTGRADUATE 
SCHOOL, Monterey, California. Use of the program at another 
installation will require appropriate modifications. 

B. ELASTIC STRESS ANALYSIS PROGRAM 

Several changes were required to add a transient thermal 
stress analysis capability to the stress analysis program 
[8], These included: (1) Addition of a consistent load 

generation subroutine LODGEN, (2) Alteration of the subrou- 
tine SOLVE to allow multiple loads, (3) Changes to subrou- 
tines INPUT, SORT, STRESS to provide steps for various 
calculations pertaining to the temperature distributions, 

(4) The capability to use nodal properties requiring the 
addition of subroutine NDPROP and numerous alterations to 
the integration process for the computation of the "stiffness 
matrix", and (5) The addition of a subroutine PRNSTR to 
calculate the magnitude of the principal stresses. 
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The program uses external storage for all information. 
Therefore, numerous additions and changes were made to the 
various WRITE and READ statements. A thorough knowledge 
of the direct-access statements used in Fortran IV is 
essential to understand the program used in this analysis. 

1 . Subroutine LODGEN 

This subroutine performs calculations used to form 

the consistent load vectors as described in Section IIA. 

\ 

A lengthy computing time is required for this subroutine due 
to the number of computations and the many direct access 
file operations. The number of Gauss points (NGPL) used 
can be specified independently of the number used for the 
formation of the "stiffness matrix". The consistent load 
vector forms the right hand side of equation (16) and must 
be of sufficient accuracy to insure valid final results. 

For the calculations made in Part I , four Gauss points were 
considered sufficient. 

2 . Subroutine SOLVE 

The basic subroutine used to solve equation (15) 
is contained in Refs. 8 and 11. In order to solve multiple 
load vectors, changes were incorporated to make necessary 
computations to all load vectors simultaneously. As des- 
cribed in Refs. 8 and 11, calculations are performed on 
each block of values in the load vector until all nodes 
have been processed. The changes which were made performed 
these calculations to a specified block for all load vectors 
prior to advancing to the next block. This method was used 
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to minimize the number of computations and the direct 
access operations. 

3 . Subroutine NDPROP 

This subroutine calculates the material properties 
at a node or Gauss point. If the elastic matrix is re- 
quested, the various properties are calculated for the 
position in question, and the matrix is formed. If the 
coefficient of thermal expansion is requested, it is calcu- 
lated for the position specified and returned to the calling 
program. This subroutine uses a separate subroutine, SHAPE, 
to evaluate the shape function values used in the calculation 
of the property values. 

4 . Preparation of Input Data 

The preparation of all data cards is explained in 
the initial comment cards of the program listing, Appendix 
G. Special care must be given to the details of the input 
cards . 

a. Object Time Formats 

The use of object time formats is made throughout 
the input subroutine. The use of object time formats elimin- 
ates the need to recall what format a certain data card is 
to follow. The object time format technique is covered in 
most books on the Fortran language. 

b. Data Check 

A data check can be made by specifying NGPS=0 . 
This requires the full 250,000 bytes of storage, but is 
completed in less than 20 seconds. 
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c. External Storage of Load Conditions 



If the load conditions of temperature distribu- 
tion and body forces are prepared by a preprocessing program, 
these values are passed using files 33 and 34. File 35 for 
use with nodal displacements was added to the input routine 
for future development. When the values are passed from a 
preprocessing program, a time value must be provided at the 
beginning of each set of load conditions. 

5 . Job Control Language 

Since the stress analysis program stores all infor- 
mation on external devices, the job control language is a 
very critical part of the program. A very thorough under- 
standing of the specific requirements for a particular com- 
puter facility is needed to use this program. All examples 
and instructions presented here apply to the W.R. CHURCH 
COMPUTER CENTER AT THE NAVAL POSTGRADUATE SCHOOL, Monterey, 
California . 

External storage operations used in the program are 
divided into two types. One type is the sequential file 
READ or WRITE statement. This method of external storage 
performs the operation on a specified record length in a 
sequential manner. The second type of operation is the 
direct access method. In this case, a specified record 
length is written into or retrieved from a given location 
within the file. In the first type of operation, no 
changes to the program language are necessary. The second 
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type of operation requires changes to the DEFINE FILE 
statements listed in the control program. These statements 
will vary according to the size of the problem being 
analyzed. 

Each file used by either a sequential type operation 
or a direct access operation must be described in the job 
control language. For the case of the direct access opera- 
tions, the DEFINE FILE statements and the job control language 
must agree both in record length and the number of records. 

In order to provide a complete summary of all files. 
Table 5 was prepared. If a particular problem does not 
require a file, its job control language and any DEFINE FILE 
statements can be eliminated. Table 6 describes each DEFINE 
FILE statement and job control card. For each file, a for- 
mula is supplied which determines the record length and the 
number of records which that file will contain. The program 
listing. Appendix G, and the job control card listing. 

Appendix D, are examples for use with a thermal stress 
problem using 278 nodes, 27 elements, and 10 load vectors. 

C. POST PROCESSING PROGRAM 

The elastic stress analysis program provides a consider- 
able amount of information. For a problem of 200 nodes 
with 10 load conditions, about 16,000 lines of output are 
generated. Interpretation of this amount of output is both 
tedious and time consuming. Usually only a few nodes are 
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TABLE 5 



FILE DEFINITIONS 
File Number - Definition 



7 


— 


Storage 


of 


Total Stiffness Matrix Blocks 


8 


- 


Storage 


of 


Element Stiffness Matrices 


9 


- 


Storage 


of 


Element Connectivity 


10 


- 


Storage 


of 


Nodal Coordinate Vectors 


11 


- 


Storage 


of 


Nodal Properties 


12 


- 


Storage 


of 


Nodal Stress Vectors 


13 


- 


Storage 


of 


Nodal Strain Vectors 


14 


- 


Storage 


of 


Concentrated Load Vectors 


15 


- 


Storage 


of 


Boundary Condition Vectors 


16 


- 


Storage 


of 


Element Material Numbers 


17 


- 


Storage 


of 


Reaction Vectors 


18 


- 


Storage 


of 


Nodal Temperatures 


19 


- 


Storage 


of 


Nodal Coordinates by Element 


20 


- 


Storage 


of 


Load Vector Blocks 


21 


- 


Storage 


of 


Nodal Displacement Vectors 


22 


- 


Storage 


of 


Predefined Boundary Displacements 


23 


- 


Storage 


of 


Nodal Properties by Element 


24 


- 


Storage 


of 


Nodal Temperatures by Element 


25 


- 


Storage 


of 


Nodal Body Forces 


26 


- 


Storage 


of 


Nodal Displacements 


27 


- 


Storage 


of 


Nodal Body Forces by Element 
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File Number 



Definition 



28 


Storage of Nodal Displacements 
by Element 


29 


Storage of Nodal Force Vector 


30 


Storage of Time Vectors 


33 


Storage of Preprocessed Temperature 
Dis tributions 


34 


- Storage of Preprocessed Body Force 
Distributions 


35 


Storage of Preprocessed Nodal 
Displacement Distributions 



TABLE 6 



DEFINE FILE Statements and Job Control Cards 
General format for these cards: 

DEFINE FILE XX (NRX ,LRWX , U , IXX) 

//GO . FTXXFOO 1 DD UNIT=SYSDA, SPACE= (LRBX, (NRX,1)) , DISP= (NEW, DELETE) 
where XX, NRX, LRBX are defined in the table. 



FILE 


Number of 
Records 


Record length 
in BYTES 


Record length 
in WORDS 


XX 


NRX 


LRBX 


LRWX 


7 


MMxNN xN RE C 7 


NCOUNTx8/NREC7 


NCOUNTX2/NREC7 


8 


4xNEL 


7200 


1800 


9 


NEL 


80 


20 


10 


NDPT 


24 


6 


11 


NDPT 


24 


6 


12 


NDPTxNLOAD 


56 


14 


13 


NDPT 


56 


14 


14 


NCLD 


32 


8 


15 


NPBC 


16 


4 


16 


NEL 


4 


1 


17 


NPBC 


32 


8 


18 


NDPTxNLOAD 


8 


2 


19 


NEL 


480 


120 


20 


NNxNLOAD 


NSx8 


NSx2 


21 


NNxNS/NDF 


24 


6 
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XX 



NRX 



LRBX 



LRWX 



22 


NPBC2 


24 




6 


23 


NEL 


480 




120 


24 


NLOADxNEL 


160 




40 


25 


NLOADxNDPT 


24 




6 


26 


NLOADxNDPT 


24 




6 


27 


NLOADxNEL 


480 




120 


28 


NLOADxNEL 


480 




120 


29 


NLOADxNDPT 


24 




6 


30 


NLOAD 


8 




2 


33 


NLOADxNDPT 


8 




2 


34 


NLOADxNDPT 


24 




6 


35 


NLOADxNDPT 


24 




6 


Notes : 










(a) For 


a 20 nodal point 


brick, 


NREC7 = 4 




(b) For 


a 32 nodal point 


brick, 


NREC7 = 6 




(c) Variables MM, NN , NS, 


etc . 


are defined 


in Ref. 8 
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considered at one time and the remainder of the information 
is reserved for later use. 

Through the use of job control language, the important 
information can be stored on external storage devices and 
processed by a separate program. A post-processing program 
was written to analyze the stored information. 

1. Program Description 

The program listed in Appendix G performs the 
following operations: (a) information stored on external 

devices is read for a specified number of nodes; (b) the 
x,y, and z components of all stress vectors are analyzed 
for the principal stress values, which are used to determine 
the maximum and minimum values for a given load condition; 

(c) Offline and Calcomp plots of stress versus coordinate 
position and temperature versus coordinate position can be 
obtained; (d) Offline and Calcomp plots can be obtained for 
the stress values in polar coordinates. 

2 . Input Data Preparation 

The input data deck has been designed to consist of 
three or four cards . The preparation of these cards is 
explained in the initial comments of the program in Appendix 
H. When more than one problem is processed, care must be 
taken to limit the number of Calcomp plots requested. All 
Calcomp plots are performed twice to make the lines dark 
enough for reproduction by other means . 
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3 . Job Control Language 



Since the information is to be passed from the 

stress analysis program, a few changes to the job control 

language discussed in Section IIIB are necessary. Files 

10, 12, 18, and 30 are changed to the following format: 

//GO. FTXXF0 01 DD UNIT=XXXXX , SPACE= (LRBX , (NRX , 1) ) , 

// DISP= (NEW, KEEP) , DSN=YYYYY, VOLUME=SER=ZZZZZ 

where XX, LRBX, NRX - defined in Section IIIB 

XXXXX - a unit number 

YYYYY - a file name, e.g. STRESS 

ZZZZZ - a device number 

The unit number and device number are obtained from the 

particular computer facility. 

The job control language for the post-processing 

program must match with that passed by the elastic stress 

analysis program. Files 8, 9, 10, and 11 are used for the 

stress values, coordinate values, temperature values, and 

time values respectively. The job control cards for the 

post-processor have the same form as those in the elastic 

stress analysis program. The linkage between the two programs 

is formed by specifying the same unit number, device number, 

and file name for the information given by the stress program 

and received by the post-processor. For example, file 12 

in the stress program would be 

/ / GO . FT12F00 1 DD UNIT=2 311 , SPACE= (LRBX , (NRX , 1) ) , 

// DISP= (NEW, KEEP) , DSN=STRESS , VOLUME=SER=SYSO 01 

and, therefore. File 8 in the post-processor would be 
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// GO . FTO 8 FO 0 1 DD UNIT=2 311 , SPACE= (LRBX, (NRX, 1) ) , 

// DISP= (OLD, KEEP) , DSN=STRESS , VOLUME=SER=S YSOO 1 

The post processing program also contains DEFINE FILE 

statements which follow the guidelines in Section IIIB. 

Further details on the procedures for job control 
language can be found in Ref. 12. 

D. SYSTEM USAGE OF ALL PROGRAMS 

The discussions thus far have pertained to the individual 
programs, however, the actual production usage of the pro- 
grams was done in a multistep process. The usage of the 
programs as a system first requires the formation of the 
system, then the production run of a problem. 

1. System Formation 

Appendix D lists the job control language used to 
create the library of programs. Details of this procedure 
are provided in Ref. 12. 

2 . Production Runs 

Appendix E lists the job control language for the 
multistep process of temperature and stress distribution 
determination. The major change in the methods is the 
transport of the temperature distribution from the first 
step to the second step. This is accomplished by the card 
beginning with //GO.FT07F001 for the first step job control 
language and the card beginning with //GO. FT33F001 of the 
second step job control language. These two cards connect 
the informatin generated and used by the two programs. 
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Reference 12 contains full details of the job control 
language used. 



IV. RECOMMENDATIONS 



As with any endeavor to construct a usable tool for 
engineering analysis, the job is never done. The following 
changes are recommended for future development: 

A. PROGRAM FOR THE TEMPERATURE DISTRIBUTION 

The present version of the temperature distribution is 
designed for an "in-core" solution technique, which restricts 
the problem size to within the computer core size. The 
stress program uses an "out of core" technique and is 
restricted only by the availability of external storage 
devices. In order to eliminate the core size restriction, 
the temperature program requires complete revision to permit 
the use of external storage devices. 

B. INITIAL STRESS DISTRIBUTION 

The load generator should be altered to allow a non-zero 
initial stress distribution. 

C. JOB CONTROL LANGUAGE 

The amount of job control language allowed at a par- 
ticular computer facility is limited. For the facility 
used during this analysis, it appears that no more job 
control language can be used with this program. The various 
files used in this program must be combined to reduce the 
number of job control cards needed. 
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D. MATERIAL PROPERTIES 



In order to more accurately solve the real engineering 
problems, the capability of temperature dependent properties 
should be added to both the temperature program and the 
stress program. 
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SAMPLE DATA DECK FOR TEMPERATURE PROGRAM 



O' 

CO 



«n 

m 

cm 

h 

o 

O' 

CO 

r^ 

o 

in 

Nf 

m 

CM 

•—4 

o 

O' 

00 

1^- 

'0 

in 

m 

CM 

r-4 

O 

O' 

CO 

vO 

m 

t 

CO 

eg 

-H 

o 

O' 

•CO 

O 

in 

m 

CNJ 

H 

o 

O' 

00 

h- 

in 

't 

m 

<NJ 

o 

O' 

CO 

h- 

vO 

in 



m 

nj 

o 

CM 

CO 

o 



OOOOOOOOoOOOOrmnSf^inNOOOOOOOOOOOOO 

Momotnoooooino^oHstN^stsooooooooooooQ 

MOMnMOOlTiOONinfNjOM'tsOM^^OOOOOOOOOOOOO 

(VO^NOinOMnocONvOiriHO'NHO'NOOOOOOOOOOOOO 

ocm sOUNUNUNuNsOuNUNvOUNuNUNUNUNsTst uNM-st mmmmmmmmmmmmm 

• CNJ 

OCMOOO OOOOOOOO OOO OOOOO OOOOO OOO OOOOO 

r- 

O' 

UN 

CO 

(NJ'O 

o 

o 

o mmr-minr-o ooo oooo ooooo 

CNJ (W |^|s t( M^N0UN0 l iN00000 i n0ir\0 

M^vOCg^OONUNMOOlAOON’NOJO 
^0'r^^Hj'r s -o < cor^'OuNoi s "UNOcor^ >0 un 

O O O O O O o O O O O O O m -h ^ (NJ p-i — 4 ro (NJ r\l c\| (\l I o c\l (NJ CO nj ni CNJ cnj 

ooooo ooooooooooooooooooooooooooo 



Cvj 

CM 

•—I'M 

(M 

cm 

CM 
I s - 
CM O' 
in 
»o 
m 
o 
oo 



oooooooo ooo 
OOOOOOOO ooo 
oooooooo ooo 
oooooooo ooo 



l/) 



ooo 

•—4 

-HC<N 

r^co 

com 

UN sO 

»-H H 

(MM- 

CMCM 

sor- 

<M(\J 

OCM 

men 

OCM 

•—4 H 

(MM* 



-^m 

CM CM 

O' r '*4 
(Mm 

r^co 



oooooooo 

oooooooo 

oooooooo 

oooooooo 



OOO O O UN UN in o Oo OOO ooooo OOOOO UN UN UN ooooo 

<m mm mmooo mm m ooooo mmmmm 



0 '- 



vor^- 



*m 



J-(M OOOOO- 
2 
UJ 

X (Nl 
UJ (M 
OM"M 
UJ CM 
CM 

O CnJ O 

h--H(Mco-tuNNO^oooo— f(Mm-tuNvor^ooo>o— KMmM’iA'Or^QO^o-Hm t 1 n nD o 



M" UN 



OCM 
CM (M 





J— 


^^^^^H^r-^r^^^cMnjcMCMcMnJcMCMCNjcMmmmnjcMCMOM’ 


m 


- | UN 


m o 


NJ 


_j m 


O' COOO 


r-4 


< m 


njm • 


o 


h- O 




O' 


LUCMO 


(\li-H «— 4 


00 


21 


m 




O 


o 


o 


UJ 


o 


UN 


l/> (M 


o 


M- 


<(Mm 


-^—4 CM UNO 


m 


com 


• UN 


CM 


•-H 


o-* 



239 



00000 0.00100 1.00000 



o 

o 

o 

o 

o 

• 



o 

o 

o 

o 

—i 

o 



«ntn m-nm uvn m 
X cm <\i <\j cnj r\j x x 
Q'OOO'OOOO 
qo n in in in in in la in oo 

sOsO^OOsO'O'OsO 
O'O'O'O'O'OON^ 
sQ vQ \Q sQ ^ >Q O sO 

A’Amininninm 

l I I l I i l I 



ininmininintnin 

oococococotocoto 

r-^ r —4 «~ ( r — 4 »— «♦ *— ♦ >H f — 4 

inininmLAinininoo 

rA'AfArA'AmcnrA 

0 'C S ^ 0 > 0'000 

00 X) rococo TO O’X) 

XxfXXfxjxXX 

oooooooon 

oooooooo I 

oooooooo 



XIXXXXXXX 

la min in m in in m 
x x rvj x x x x x 
OC 

in in in in in n m m 
xxxjxxxxx 
OO'C^OOOOO 
mminminminm 

xxx XI x X X x 

I I I I I I I I 



O 

40 



■ n in in in in in in in o 
oo oo o ' oco oo oo oo o 

i—4 rH iH r -4 fH »-“4 O 

co m n innn in m in 
n m in <n a acapa 
o > t> xc v cnc^ 
rococo o to to co co 
X X X XX X! X x 



in in in in in in in in 

X) X X XI X X X X 

cm> o> > o* c* o o* 

• uo la n m uA in in n n co 

-4 

O 0 *0 OOsO o o 
0 > 0 ' 0 S 0 ' 00 ' 0 ' 0 ' 
O'O 0 * 0 v 0 v 0 o o 



^OOOOOOOO'tr 0 Nl s, r s, l s, NN^I s » oA 
I OOOOOOOO •••••••• 

• » • • • • * • n in in in m uvn m 

oooooooo i i i i i i i i c 



(xj-HvOCM^-coinostPATOinOvt- 

r-<XXX~4 (XI XI XI r-^ 



cxj xia- vt 'OincorAO'Axjr^vt omco rAO r 

H l'AX|AJr “"4 *“H fAXIXIr -4 



4 (\J *■ 



o CM^ co in O Nt- ' A CO 

•— H XJ XJXI i —4 XJ 



oc^r^ooinovr^ooino-t— »'O r >r^xix;r^<rsoinco(Aa>xixii^^ omcorAOrA 

~*XJ (M(AJr-l (XlXIXI—f —4 'A XIX ~4 •— < pT) XX ~4 



-H XXXr-W XI 



240 



0 Q 2 li 9 jt»id 51 



in m m mm m m 
CO CO 00 00 CO 00 00 

mmmmmmm 
canmncn nm 
0 s ( 7^0 O'Q'O o 
CO 00 CO CO CO 00 CO 
f\JC\]C\]r\lC\jnirNj 
ooooooo 
ooooooo 
• •••»•• 
ooooooo 



mo>t * 

f\JC\|r-< 



r\jrvjcsjr\»(Mf\jr\i(\j mminmininmin 
o a> cn o o o o o ooaoooooaocococo 

in nmmmmmm 

go'xi m nj njrvjnj m eg com mm nmmmm 
cn n ^ o ^ cn 

m m m m m minn a* m c* m o o ct* O' 
(\injc\jr\t<\jc\j(\jc\j co ooo oco ocqco 
nnOOMXiOO'Cri (Sjc\j<\|r\|njnjc\jr\j 

rommmmmmmmroooooo ooo 

I • I OQOOOOOO 

C\|C\l<\]iXJr\j<\J AJCvJ 

t t 1 I I t l t oooooooo 
4 o eh <m c\j n- so m oo m avn v4» go m co n cr» 

—4 r00sJC\J-H *-H ♦ mOsJCM^H r-H 



movr^-ocnr^njcMr^sf om coin cnc\jf\n^<r ^Qm corner 
corn-* ^ m c\jr\j— ♦ «-h *— i rn ( Njnjr-i ^hhci 

Q 

b- 

co 



241 



! 







O' 












CO 












r- 












o 












LA 
























cq 












<M 












•—4 












o 












O' 






oo 






00 






00 






I s *. 






< 






NO 






Q. 






in 






•* 












3C 






rq 






UJ 






eg 






Z-*-* 






•^4 






— LULU 






o 






II h-H- 






O' 




Q 


CL UJlU 






00 




QC 


10-J-J 






h- 




< 


— UJLU 






v0 




o 


QQQ 






in 
















■* 


Q-3: 3 






cq 




\ 


s:ujuj 






eg 






UJ2!Z 






rM 




> 






O 


o 




CO 


O i| || 




zs 


O' 






II CL CL 






00 




Q 


Zoo 00 




ho: 






UJ 


L/) — — 




ooo 


o 






QQQ 




hO 


in 




O 


•* •* 


CO 


-JCZ 






o 






a. 


cq 




_J 




X 


Q 


c\) 




o 


HHH 


— 


aC lu 


'■'■4 




ll 




Q 


< aC 


o 






CMrgcM 




03 


O' 




z 




Ui 


h- 


CO 




< 




CL 


o< 


h*i 




a: 


j^i j 


Q- 


Ocz 


v0 




o 


>>> 


< 


OC UJ 


in 




o 


ooo 




l-CL 






a: 






Z2 


rq 




a 


u ir ii 




Qlu 


CM 






LULULU 




Oh 


•—4 




UJ 


ooo 






o 




a: 


<<< 




cQa: 


O' 




3 


CL CL CL 




oo 


CO 


O 


h- 


oooolo 




">LL 


N. 


n 


< 








v0 


cq 


OS 


<x<< 






in 


II 


UJ 


OOQ 








o 


CL 


LOOOOO 






rq 


o 


Z 


>>> 






CM 


• 


LU 


oooooo 






r-4 


z 


h— 


li ii it 






O 


o 










O' 




LL 


— — — 






00 


o 


O 


zzz 






r-- 


LU 




333 






O 




UJ 








n 


O - 


Q 


OQQ 






M- 


oCOQ 


O 


QQQ^ 






cq 


<oo 


o 








CM 


oo 




HfHHC 



o 

O' 

CO 

h- 

L A 
<t 

cq 

00 




Q 

CC 

< 

O 



> 

CO 

Q 

LU 

!£ 

o 



3 

LL 

X 

O 

UJ 

Q 



< 

o 



242 



SAMPLE DATA DECK FOR STRESS PROGRAM 
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S SUBROUTINE IS THE BLGCK SOLVER PATTERNED AFTER A ROUTINE 
GINALLV WRITTEN BY PROF. J.G. CANTIN AT THE U.S. NAVAL 

TGRADUATE SCHOOL 
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SUBROUTINE PRNSTR ( STR , Si t S2 ,S3 ) PSTOOOIO 

THIS SUBROUTINE DETERMINES PRINCIPAL STRESS VALUES BY CUBIC 
EQUATION FORMULA 
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